Spatio-temporal distribution of nucleation events during crystal growth 
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We consider irreversible second-layer nucleation that occurs when two adatoms on a terrace meet. 
We solve the problem analytically in one dimension for zero and infinite step-edge barriers, and 
numerically for any value of the barriers in one and two dimensions. For large barriers, the spatial 
distribution of nucleation events strongly differs from p^, where p is the stationary adatom density 
in the presence of a constant flux. The probability Q{t) that nucleation occurs at time t after the 
deposition of the second adatom, decays for short time as a power law [Q{t) ~ t"^/^] in d = 1 and 
logarithmically [Q{t) ~ l/ln{t/to)] in d = 2; for long time it decays exponentially. Theories of the 
nucleation rate u) based on the assumption that it is proportional to are shown to overestimate 
a; by a factor proportional to the number of times an adatom diffusing on the terrace visits an 
already visited lattice site. 



Molecular Beam Epitaxy (MBE) is one of the most 
common techniques for growing nanoscale materials [Q: 
a controlled flux of particles arrives ballistically on a sub- 
strate whose temperature is generally high enough to ac- 
tivate surface diffusion of newly deposited adatoms. fn 
the absence of preexisting steps, i.e. for high-symmetry 
surfaces, growth proceeds through the formation of stable 
dimers (or largest nuclei, depending upon the tempera- 
ture, the flux and surface symmetry [^) and the subse- 
quent aggregation of diffusing particles. 

Nucleation processes take place on a flat surface during 
the submonolayer regime Q and afterwards on islands, 
fn this Letter we are interested in the latter case, where 
the process of nucleation occurs as follows: adatoms are 
deposited at a rate F per unit time and lattice site and 
diffuse at a rate D. The typical time between two depo- 
sition events on a terrace of size L is Tdep = {FL'^)~^ and 
the typical residence time of an adatom on the terrace 
is jj] Tres ^ L[L + OLdt-Es) I F) , whcrc ^es is the Ehrlich- 
Schwoebel (ES) length % that measures the strength of 
step-edge barriers 1^ hindering the descent of adatoms; 
Old is a numerical factor depending on dimension only. 
A nucleation event takes place when a newly deposited 
adatom finds the previous one still on the terrace and 
they meet before descending, forming a stable dimer. 
Once two adatoms are on the terrace at the same time 
a nucleation event takes place, if it does, in a time of 
order Tt,- ^ jD, the typical time required for a traver- 
sal of the terrace Q|. Therefore the probability that a 
third atom influences the nucleation process is of order 
P3 ~ Ttr/Tdep^ which is vcry small in all physically rele- 
vant cases. More precisely, in d = 2 ~ L'^ / where 
£d ~ [D/Ffl^ - the so called diffusion length |] - is the 
typical size of a terrace when nucleation occurs. Since 
L < £o and 3> f we conclude that ^ 1 and for irre- 
versible nucleation only processes involving two adatoms 
must be take into account. 



The simplest theoretical treatment of nucleation phe- 
nomena in crystal growth is based on the assumption 0| 
that the nucleation rate to is proportional to the square 
of the stationary adatom density p in the presence of a 
constant flux F. Such a mean-field treatment has been 
used to estimate physical parameters of materials 
and as an ingredient in mesoscopic models of epitaxial 
growth However, recent papers [^,0 have shown 

(via Monte Carlo simulations and scaling arguments ll^] 
or theoretically that the mean field prediction uj ~ 
for the nucleation rate is not correct in d = 2 for large 
ES barriers. This prompted us to reconsider the general 
validity of mean-field theory, ft is a priori not clear to 
what extent two-particle properties as uj can be simply 
obtained from p, a quantity describing a single particle on 
a terrace. In this Letter we solve exactly the problem of 
irreversible nucleation on a terrace in I and 2 dimensions 
for all values of the ES barrier. We compute (analytically 
or numerically) the spatial and temporal distributions of 
nucleation events and the total nucleation rate. In this 
way we are able to assess when and why the mean-field 
approximation provides sufficiently accurate results and 
by how much it fails. 

Let us start with the case ^es = 0. The discrete evolu- 
tion equation for the probability Pn{t) of finding a single 
adatom in site n at time t is 

Pn{t+l) = ^[pn+l{t) + pn-l{t)]. (I) 

The boundary conditions are po(t) = pL^i{t) = 0. 

By looking for solutions of Eq. (|^) of the form p„ (t) = 
N„T{t) one finds the general solution 



p„(i) = ^ Afc cos* (fcn) sin {nkU) , 



(2) 



fc=i 



with coefficients = Pn(0) ^iii('^^n), where 

U = Tr/iL + l). 
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Nucleation occurs when two adatoms wandering on the 
same terrace meet. The one-dimensional diffusion of two 
particles can be reformulated by taking the coordinates 
TO, n of the two walkers as the coordinates of one walker 
in d' = 2 moving on a square lattice {L x L). The gen- 
eralization of to two dimensions is Pm.n{t + 1) = 

j[Pm+l,n{t) + Pm~l,n{t) + Pm.,n+l{t) + P,nji^l[t)] with 

boundary conditions po,n = PL+i,n = Pmfi = Pm,L+i = 
0. Again, separating space and time variables one obtains 
the most general solution 

Pm,n{t) ^ J2 ^ [Cos(fcn) -fCOs(in)]* 

X sin(mfcn) sin(nj'n), (3) 
where the coefficients Bkj are 

■^''^^ (lTt) ^ Pm,n(0)sin(TOffl)sin(njn) . (4) 

The initial condition Pm,n(0) = PmPn encodes the in- 
formation about the way particles are deposited. Notice 
that t = indicates in Eqs. the time when the sec- 
ond adatom arrives. We consider adatoms landing on 
the terrace with spatially uniform probability = 1/L. 
However, since arrivals on the terrace are not simulta- 
neous, the distribution of the first adatom has changed 
to Pm{t') (Eq. H) when the second one lands on the ter- 
race. The nucleation process depends therefore on the 
precise interarrival time t' , which is a poissonian random 
variable ParhX^') = [rdep)'^ exp{-t' / Tdep) ■ 

One must compute physical quantities O by averaging 
over t': O = J^^^o PAnn{t')0{t'), where 0(V) is com- 
puted using as initial distributions Pm{t') for the first 
particle and p^ = 1/L for the second. If 0{t') is lin- 
ear with respect to the initial distributions (as are all the 
quantities considered below) , the sum can be replaced by 
a single computation of O with the initial distribution of 
the first adatom given by pm = Z]t^=n_fARR (i')Pm(^')- 

The general solution can be found |l2] , but in the limit 
Tres ^ Tdep, which is rcalistic for MBE, one can consider 
-PARR(i') as a constant on the scale Tres and then pm — 
'^t^=oP"iit')- Therefore is proportional to the nor- 
malized stationary solution p^ of (Q) in the presence of a 
constant fiux p^ = 6/[L{L + l){L + 2)]m{L + l-m)] 
will be taken as the distribution for the first particle 
for computing P{n) and Q{t). 

So far we have not considered the interaction between 
adatoms, i.e. the fact that when the two particles meet 
they stop diffusing. An irreversible nucleation event gen- 
erally occurs when two adatoms are on nearest neighbor 
sites. To allow an analytic treatment we consider here 
a nucleation event to occur when the adatoms are on 
the same site. When this occurs the dynamics stops, 
implying Pm.m 

(t) = for all t > 0. This boundary 
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FIG. 1. Main: normalized distribution of the nucleation 
sites P{n) in d = 1 for ^bs = and L — 100. pn ~ is the 
(normalized) stationary solution of Eq. (|l|) in the presence of 
a constant flux. Inset: The same for ^es ~ oo. 

condition can be taken into account by the classical 
image method |p^ , |l5t : The initial condition is chosen 
to be antisymmetric with respect to particles exchange 
Pm.niO) = -Pn,m{0)- This implies Bkj = -Bjk and 
therefore Pm.nit) — —Pn.m{i)- In this way the bound- 
ary condition Pm,m{t) = is obeyed for all t because 
the two triangles (to > n) and (to < n) are dynamically 
disconnected. 

When the antisymmetric initial condition is not im- 
posed adatoms diffuse without interacting: they do not 
feel each other even if they are on the same site and 
they wander until they get off the terrace. In the follow- 
ing we will speak of nucleations for these noninteracting 
adatoms, intending that they "nucleate" when they are 
on the same site. Clearly two noninteracting adatoms 
may give rise to several "nucleation" events before leav- 
ing the terrace. 

For noninteracting adatoms the coefficients B^j are 
simply the product of the single adatom coefficients 
AkAj. In the case of interacting adatoms we ob- 
tain (l| Bk, = [2/{L+l)f [B< - B<], where B< = 
Em<n-Pm:n(0) sin(TOm) sin(njn). 

For interacting adatoms the probability of a nucleation 
event on site n at time < + 1 is given by (l/2)[pn,n+i{t) + 
Pn-i,n{t)], while for i = it is Pn,n{0). For noninteract- 
ing adatoms it is simply Pn,nit) at any time. By sum- 
ming over t one obtains the spatial distribution P{n) of 
nucleation sites, which is reported in Fig. |l|. It turns out 
immediately that the mean-field prediction is exact 
only if adatoms do not interact, i.e. nucleation events 
following the first one are taken into account. In the in- 
teracting case the distribution of nucleation sites differs 
from but the discrepancy is rather small. 

We now discuss the distribution of nucleation times, 
i.e. the probability Q{t) that adatoms meet at time t 
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FIG. 2. Main: Q{t) on a log-log scale for L = 100 in d = 1 
for l-Bs ~ (bottom) and £es ~ oo (top). The straight line 
shows a decay of Q{t) as t~^^^. Inset: 1/Q{t) on a log-lin scale 
for L = 40 in d = 2 for ^es ~ (top) and ^es ~ oo (bottom). 
The straight line shows a decay of Q{t) as l/ln(t/to). 

after the deposition of the second adatom. 
For non-interacting adatoms, 

^ L + 1 ^ 

Q{t) = ^p„,„(i) = ^ Bkk cos*(fcn) . (5) 

n=l fc=l 

We can rewrite cos*(0) — exp[t lncos(0)] and since 
the coefficients Bkk diverge for small k as we 
expand the cosine for small and finally extract 
the dominant contribution coming from the mode 
k = 1: Q{t) ^ exp[— In the case of in- 
teracting adatoms Q{t + 1) = l[Pn,n+i(*) + 
Pn-i.nit)]- By evaluating the summation, we ob- 
tain Q{t + 1) = Eij^i BkjCkj2-' [cos(m) + cos(jn)]*, 

where the Ckj = X)n=i sin[fcn(n — 1)] sin(jnn) can be 
evaluated exphcitly. Approximations along the lines of 
the above treatment lead to 




The sum can be rewritten as the integral exp[- • ■]dk 
plus the boundary term for k — \: (1/2) exp[— (5/4)n'^t] 
(the boundary term for k = L is always negligible). It is 
easy to see that the integral prevails for t <C L^/tt and 
gives Q{t) ~ 16\/2/ {LT:^^^\/t), while in the opposite limit 
we have an exponential decay. The existence of the two 
regimes is clearly shown in Fig. ^. The two behaviors 
can be interpreted physically. For short times, terrace 
edges can be neglected and one can focus on the relative 
coordinate (n—m) of the two particles. Nucleation occurs 
when (n — m) vanishes for the first time and Q{t) ^ 
t~^^'^ is simply the spatial integral of the first-passage 
distribution probability u&. The exponential decay for 
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FIG. 3. Main: values of the normalized distribution of nu- 
cleation sites along the diagonal P{n, n) for Ies = 0, d = 2 
and L = 20. Inset: the same for ^es = oo. 

long time is the effect of the decaying probability that 
the adatoms remain on the terrace times the vanishing 
probability that they have not yet met. 

In two dimensions the two diffusing adatoms can be 
mapped into a four dimensional problem for a single ran- 
dom walker: Pmi.ni.m2,n2i't) is the probability of finding 
one atom on site {mi,ni) and the other in (m2,?T.2) at 
time t. The solution in the noninteracting case on the 
four dimensional hypcrcube is easily found. However, 
the simple generalization of the interacting case to d = 2 
is not possible because the plane corresponding to nu- 
cleation events does not divide the hypcrcube into two 
dynamically separated regions, reflecting the fact that in 
two dimensions two adatoms can exchange their positions 
without meeting. The results shown for the interacting 
case are obtained by numerically solving the evolution 
equation for the probability p. 

The main part of Fig. ^ clearly shows how the re- 
sults obtained in d = 1 for the spatial distribution of 
nucleation sites are also true in two dimensions: for 
non-interacting atoms P(m,n) = n'j the inter- 
acting case a discrepancy exists but is practically neg- 
ligible. The distribution of nucleation times Q{t) de- 
cays exponentially for long time, as in d = 1. Again, 
at short time it can be derived ||l^ by using the first- 
passage probability arguments, yielding -for interacting 
adatoms- Q{t) ~ l/lii{t/to) (Fig. |, inset). 

Let us now consider €es = oo. In this case, the 
boundary conditions for a single adatom in d = 1 
are poit) = pi{t) and PL+i{t) = PL{t). The gen- 
eral solution is now a superposition of the functions 
Xk{n) — tan[fc7r/(2L)] sin(nfc7r/L) -I- cos(rifc7r/L), with 
fc = 0, . . . , L — 1 and the solution of the nucleation prob- 
lem is found much in the same way as for zero barri- 
ers ]T^ . The initial condition is now with both adatoms 
uniformly distributed, because for ^es = oo the station- 
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ary solution for the single adatom is = 1/L. 

The spatial distribution P{n) of nucleation events is 
compared to the mean-field prediction {P{n) = l/L'^) in 
the insets of Fig. |^ (d = 1) and Fig. || (d = 2): they are 
in striking disagreement. We conclude that -concerning 
P{n)- mean field theory is a good approximation for 
weak step-edge barriers, but the agreement is extremely 
poor for strong ones. 

The temporal distribution Q{t) of nucleation events 
does not change qualitatively when the value of t 

ES IS 

varied (see Fig. 2). The power-law {d = 1) and logarith- 
mic [d — 2) decays at short time do not depend on the 
finite size of the terrace and therefore on £es j which enters 
in boundary conditions only. The typical time necessary 
for the adatom to feel the presence of the boundaries 
corresponds to the time Ttr ^ l? jD. For longer times 
Q{t") vanishes exponentially. In this case the probability 
for adatoms to remain on the terrace is clearly constant. 
Nonetheless Q{t") must decay rapidly since its integral is 
the total probability of a nucleation event, which is 1 
for infinite barriers: The exponential decay is due to the 
vanishing probability that the two adatoms have not yet 
met before. 

^From the experimental point of view, a relevant quan- 
tity is the nucleation rate tu, defined as the number of 
nucleation events per unit time on the whole terrace. 
Within our framework, we can rigorously show, in all di- 
mensions and for all values of the barriers, that the mean 
field value Wmf = DL'^p^, with p = Fxres the average 
density of occupied sites, must be corrected by a factor 
proportional to AT = where N^^^ all are, respectively, 
the number of distinct or all sites visited by a single atom 
on a terrace. For reasons of space we only sketch here 
the derivation of this result, that will be presented in de- 
tail in a longer publication [Q. The nucleation rate lo 
is equal to the deposition rate of a single particle 1 / Tdep 
times the probability Pnuc that such a particle nucleates 
a dimer before leaving the terrace. Pnuc is the average 
over all interarrival times t' of the probability that two 
particles meet, when the second is deposited a time t' 
after the first. As discussed previously for the generic 
quantity O, Pnuc can be computed by considering the 
two particles deposited simultaneously, but with the first 
distributed spatially as Pm- This last function is shown 
to be equal to {jres/Tdep) Pm^ whcre p^ is the normalized 
solution of the stationary diffusion equation for a single 
particle. Hence we have Pnuc = (Tres/Tdep) W, where 
W = X)t^o probability that two adatoms, de- 

posited at the same time with normalized distributions, 
meet before leaving the terrace; W can be shown to be 
proportional to N^is/L'^- Putting all together and using 
Tdep = 1/{FL'^) and Fr^es = P, we obtain uj ^ FL'^pN^is- 
In a perfectly analogous way one can show, for noninter- 
acting adatoms, cjni ^ FL'^pN^ii- By considering that 
for noninteracting particles Tres = Ni^n/D one sees that 



cJni is proportional to uj-mf, so that finally 

This result holds in all dimensions and for any barrier. 

The value of the correction factor Af depends of course 
on d, L and £es- The numerator iVaii is just proportional 
to the average density of adatoms on the terrace: iVaii ~ 
L{L + a^^Es)- The value of the denominator iVdis is well 
known [l^ in absence of step-edge barriers, being of order 
L in d = 1 and of order L^/ In L in d = 2, and it is trivial 
in the limit of infinite barriers, being exactly equal to L'^. 
Hence in d = 1 we obtain Af ^ {L + ai^Es)- In d = 2 it is 
possible to find [ pTf an interpolation between the limits 
Af ^InL (4s < L) and Af - 4s/i (4s > i): - 
(l+a24s/i)M, where - l-[l-l/(lnL)]-'-/-"-. The 
previous interpolation is in reasonable agreement with 
Monte Carlo simulations |l2j performed for intermediate 
barriers. 

In conclusion we have provided exact results for the 
spatio-temporal distributions and total nucleation rates 
for irreversible nucleation on terraces during crystal 
growth. These results should be used wherever mean- 
field approximate results were so far commonly used. 
One example is the experimental determination of ES 
barriers from the rate of nucleation on terraces [^,^. 
Another example is the modelization of epitaxial growth 
with mesoscopic models : the location of new terraces 
must be chosen according to the correct spatial distribu- 
tion of nucleation events, derived here. 
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